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Abstract : 

In Part I of this paper a new hybrid-stress finite element algorithm, 
suitable for analyses of large quasi-static deformations of inelastic solids, 
is presented. Principal variables in the formulation are the nominal stress- 
rate and spin. As such, a consistent reformulation of the constitutive equa- 
tion is necessary, and is discussed. The finite element equations give rise 
to an Initial value problem. Time integration has been accomplished by Euler 
and Runge-Kutta schemes and the superior accuracy of the higher order schemes 
is noted. In the course of integration of stress in time, it has been demon- 
strated that classical schemes such as Euler's and Runge-Kutta may lead to 
strong frame-dependence. As a remedy, modified integration schemes areproposed 
and the potential of the new schemes for suppressing frame dependence of 
numerically integrated stress is demonstrated. The applicability of explicit 
and implicit forward gradient schemes to improve the stability of time- 
integration in large deformation problems is investigated. The feasibility 
and performance of the present methods are demonstrated in a number of problems, 
and it is found that the stresses obtained by the present method are of excep- 
tional accuracy, much more than could be expected of an assumed-velocity based 
finite element algorithm. 

In Part II of this paper, the topic of the development of valid creep 
fracture criteria is addressed. Until now, the so-called C* integral, intro- 
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duced by Goldman and Hutchinson, was attempted to be used in the literature 
to correlate creep crack growth rate. In the present work, a new path-in- 
dependent integral parameter (^) which has a considerably more degree of 
generality and validity than the C* integral, is introduced. The mathematical 
aspects of this parameter are first reviewed by deriving generalized vector 
forms of the parameters (T) and using conservation laws which are valid for 
arbitrary, three-dimensional cracked bodies with (macro) -crack surface tractions, 
body forces, inertial effects and large deformations. Two principal conclusions 
are that (^) is a valid crack-tip parameter during nonsteady as well as steady- 
state creep and that (^) has an energy rate interpretation whereas does not. 
Using this new integral, and a finite element analysis procedure, several funda- 
mental aspects of creep crack-growth are studied numerically. Specifically 
numerical results are presented for a double-edge-crack specimen for which ex- 
perimental results are available. Finally, a simplified methodology for pre- 
dicting creep growth behaviour is presented, based on the conclusions drawn from 
the present numerical simulations of experimental data. 

Part I ; Stress Analysis of Inelastic Solids Using Assumed Stress Finite 
Elements 

Nomenclature ; 

C(t) configuration (image of the body in space at time t) 

X position vector in space at time x 

jc position vector in space at (present) time t 

V = e^^ o/3x^; ^ 3/3X^; a = 3a/3t + V.Va material derivative of 'a* 

^ deformation function; maps C(x) to C(t) as 

X = 

V velocity function; related to deformation function as 

v(x^(x,t),t) = i^(x,t) 
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T 

H (V^^(X, t)) deformation gradient 
J = det F 

X -T 

T 

L = (V^(x, t)) velocity gradient 
j = V.^(3c, t) dilitation 
1 T 

e = ^L+L ) stretching 
1 T 

0 ) = -^(L-L ) spin 

T true traction: T nominal traction relative to C(x) 

— 

X true stress: o^= J x Kirchhoff stress relative to C(x) 

’ ~X X- 

t = F 0 nominal stress relative to C(x) 

~x ~x -X 

T true traction rate; T nominal traction rate 
a = jx + X Kirchhoff stress rate 
t = -(e+u) . X + a nominal stress rate 
0 *=g_y,x+x.u ' CO rotational ' stress rate 

Introduction ; 

The research which produced the present hybrid stress fintie element 
algorithm was motivated by the observation that hybrid stress algorithms con- 
sistently outperform those using velocity (or displacement) as the sole 
variable. Hybrid stress models for infinitesimal deformation of shells and 
incompressible solids have been topics of intense research since Plan's first 
presentation of such a model In [1]. However, hybrid stress models for finite 
deformations have only been researched since de Veubeke's [2] presentation of 
a complementary energy principle for finite elastic deformations, and Atluri's 
[3], [4] generalization of that principle for inelastic solids. A hybrid stress 
model for finite elastic deformation was presented by Murakawa [5]. In this 
report a hybrid stress model for finite Inelastic deformation is detailed. 
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T he Boundary Value Problem ; 
Compatibility 

VX(e-m) = 0; 

Linear Momentum Balance (LMB) 

V.t + = ^; 

Angular Momentum Balance (AMB) 

• • n 

[ ( e+m) . T + t ] - [ ( e+m) . t + t ] 


Constitutive Equations 




e = 9R/3r; r = — (t+T.aj-m.x+t ); 

Velocity Boundary Condition (VBC) 

6^. (-e+o)+W) = 0 on (6s^ is any tangent on S^) 

Traction Boundary Condition (TBC) 
n . t = T on S 


(I.l) 


( 1 . 2 ) 


(1.3) 


(1.4) 


(1.5) 


( 1 . 6 ) 


- ~ -t CT 

Above are listed the equations of the general boundary value problem 
associated with quasistatic deformations of inelastic solids. From (I.l) through 


(1.4) one may obtain 18 scalar equations for the 9 unknown stress rate com- 
• i-i 

ponents t , 3 unknown spin components m , and 6 unknown stretching components 
. It is possible to reset (I.l) through (1.6) so that only velocity com- 
ponents appear as variables. Alternatively one may use (1.4) to eliminate 
e as a variable in (I.l) and (1.3), thus obtaining a boundary value problem 
for the components of stress rate and spin. Any solution of this latter boundary 
value problem is a stationary point of the functional 


TT (V,m,t) = / {-R - T:(m.co) + t:m}dV 

me - - - 2 


r n , t . VdS + / (n.t-T ) 


.VdS 


(1.7) 
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( 1 . 8 ) 


provided that only stress rate variations 6t such that 
V . 6t = i) 

and spin variations such that 

6to + = 0 (1.9) 

are admitted to The velocity in plays the role of a Lagrange mul- 

tiplier. 

We replace the boimdary value problem for t and u by the generalized b.v.p. 

6tt (V,o),t) = 0 (1. 10) 

me — - - 

with subsidiary conditions: 


0 ) + 0 ) =0 
V . t + = 0^ 


(1. 11) 

( 1 . 12 ) * 


The stationary conditions of it , when (1. 11) and (1.12) decide admissible 

me 

y and t, are 


L 

L 

Is 


{ [(e+y) . T + t] :6y} 
{ [-e+uj] : 6t}dV + 


(n.t-T ) . 6VdS = 0 

— f* ““ 



a 

where e is written for 9R/3r. 
tional statements of AMB (1.3) 
(1.6), respectively. 


= 0 


n. dt.VdS = 0 


Equations (1.13), (1.14), and (1.15) 
, compatibility (I.l) and VBC (1.5), 


(1.13) 


(1.14) 


(1.15) 


are varia- 
and TBC 


* the general solution of this equation is known as t=t +c , where t°=VX0 and 
V.t^=-pb. 
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Reformulation of the Constitutive Equation : 

As may be surmised from the development above, the present approach 
necessitates reformulation of the constitutive equation. In applications 
one typically Is given or may find a constitutive equation of the form 

a* = V:e + E (1.16) 

where V and Z generally depend upon e, t, and scalar Invariants*. If V and 
Z do not depend on e then we say (1.16) Is affine with respect to e. This 
Includes most material models found In the engineering literature. If (1.16) 
Is also an Isotropic function then V and Z may be set In the forms 

'yki " 

X^^(s 6 ) + X^^(s t' ) + X^^(s s ) 

Ij kl Ij kl Ij kl 


12 3 

and Z..=n6 +nT'+vs.. 

ij Ij Id 


(1.18) 


where t' Is the stress deviator and s=t'.t'. The and are functions 

of the scalar Invariants of the deformation. 

This form still Includes most material models found In the engineering 
literature. Notable exceptions are models for materials with anisotropy In the 
stress free state, where the present model reduces to 

g* = 2y^e + A^^(l:e)I (1.19) 

By retaining the constitutive equation In the general form (1.16), (1.17), (1.18) 
In our development we are led naturally to a 'xinlfled numerical procedure’ for 
problems of large strain elasticity, plasticity, vlscoplastlclty , and creep. 


temperature, strain histories, as well as joint Invariants of J and e. 


310 



The desired form for the constitutive equation may be found by noting 
the relation between the stress rates o* and f; 


r = a* - y(t. e+e. x) 


( 1 . 20 ) 


Using (1.16) to eliminate d* from (1.20) gives 


r - (V: e+r) - — (x.e+e.x) ■= W:e + E 


( 1 . 21 ) 


where W is defined by 

W..,, = V,,,, - i (x., 6, X, .) 
ijkl ijkl 2 xk Ij Ik Ij 

Inversion of the relation (1.21) yields 


e = W ^:(r-E) 


( 1 . 22 ) 


(1.23) 


If W ^ is symmetric (i.e., then a potential R exists for e: 

e = 3R/3r; R = |-(r-E) ;W"^: (f-E) (I.2A) 

The condition necessary and sufficient for W ^ to be symmetric is 
This condition is satisfied by most engineering materials*. 

In practice one must construct W from V, then invert W (if possible) to 
achieve the form (1.24). This is a major undertaking from a computational 
point of view since W is generally different at each point of a stressed 
body. Therefore special attention is given to practical methods for construc- 
tion of W For plane problems W ^ can be found analytically. For general 
problems in which V ^ is known a simple approximation for W ^ is often of 
acceptable accuracy. The details of the two special cases are discussed in 
Appendix B of [6]. 

The Finite Element Algorithm ; 

Equations (1.13), (1.14), and (1.15) are the basis for the finite element 
algorithm presented in this report. The finite element equations are obtained 


* Flow laws using the corotatlonal rate of the true stress are an exception. 
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by Introduction of polynomial representations for V, and t* to (13) and 


,th 


(14); on the N element: 

N : isoparametric shape functions 


NQ j 

V = r N. % 
i=l 


NW 


oj = 1 QW.ot^ where QW.+QW, = 0 
~ i=l -i N ~i ~i 


NT , 

t = 1 QT.3^ + t where QT.=Vx4’ 

i=l ^ 


1 -i 

b 


(1.25) 

(1.26) 
(1.27) 


V.t“ = -Pb 

The representations for m and t are independent on each element, so (1.15) must 
be replaced by the *interelement traction reciprocity' relation; 

NELM 


NELM [ - . /• ^ ) 

2 I / (n.t.6v)dS -I T .6vds[ = 0 

N-l ) 


(1.28) 


which includes (1.15). The finite element equations are obtained by performing 
the assigned integrations (Gaussian quadrature rules are used) . Those equa- 
tions are listed below (the element index 'N' has been suppressed on the 

N 

(1.29) 

(1.30) 


and 6„) : 
N 


P 

P 

+ [G]{qj^}} = 0 


NELM mm 

r oX> - “ 

N=1 


(1.31) 


* Mathematical 'rank' conditions require that NT^NQ-T, where T is the number of 
translational degrees of freedom of an element. Moreover, the [QW] and [QT] 
should be of the same polynomial degree. 
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Henceforth we refer to (1.29) as AMB, to (1.30) as compatibility, and to 
(1.31) as TBC. The individual matrices are defined below: 


- 

r {(T.QW^):p:(x.QW.) + x:(QW^ QWj))dV 

((.32) 


J* Ux.QW2:p:(QTj) - QW^tQT^ldV 

(1.33) 

II 

1" {(QT^):D:(x.QW^) - QT^;QW^}dV 

(1.34) 

22 

f {(QT ):D:(QT.)}dV 
/ V - X = - J 

N 

(1.35) 

■ . 

f n.(QT ).(N.)dS 

Js„ ^ 

(1.36) 


T .(IL)dS 
~t — i 


.e.b ^ r 


{(QT^):(-D:t°) }dV 


’i’^ = J {(x.QW^) :D;E}dV 


= r 


{(QT^):D:pdV 


and p is obtained from W by symmetrization: 
°ljkl ^ '^jlkl ^ijlk ^jilk^ 


(1.37) 

(1.38) 

(1.39) 

(1.40) 

(1.41) 

(1.42) 


*The last term in the integrand is a residual whose significance is explained 
in [6] . 
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This symmetrization is easily done after W is computed, and serves to re- 
duce by a factor of four the number of multiplications required to compute the 
H matrices (1.32) through (1.35) , and other matrices involving W 

The procedure which leads one from equations (1.29) and (1.30) to the 
element stiffness matrix is virtually identical to that of Pian [1]. We de- 
fine the element ’H-matrix* as 


[H] = 


r-Rll jjl2-, 


21 22 
H 


(1.43) 


b Z 

and loads {P } and {P }, due to body force and fluidity, respectively, as 


{P^ = 


,a,b 


,6,b 


: {P } = 






(1.44) 


Then (1.29) and (1.30) may be collected into a single equation as 

If D is symmetric, that is, if » then from (1.32) through (1.35) we 

easily determine that [H] is symmetric. 

If the H-matrix is not singular, then we solve the matrix equation (NQ+1 
right hand sides) 



0 

G 


{q} + {P^ + P^} 


(1.45) 


[h][h“^g 


H“S] 



0 




=■ 


P + P 


G 



(1.46) 


on each element. Explicit calculation of the inverse of [H] is not only un- 
necessary, but substantially increases the time required to generate the ele- 
ment stiffness matrix. According to (1.45), the spin and stress parameters on 
each element are given by 


“1 = [H ^G]{q} + {H"^P} (1.47) 

P 
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Using (1.47) to eliminate {a/g} from TBC (1.31) leads to 

NF.LM T -1 T 

^ + {6qj^}[0 G^]{H P^} - {dq^^} {F^}} = 0 (1.48) 

N=1 

in which the element stiffness matrix has been identified as 


[^] = [0 G^][H-^Gjj] (1.49) 

and the resultant nodal ’forces’ are given by 

-[0 gL{H“Vj + {FJ (1.50) 

N N N 

It is easily verified that the element stiffness matrix [K] is symmetric if 
[H] is, and so the symmetry of [K] ultimately depends upon the symmetry of the 
constitutive matrix W. 

To this point all of the finite element equations are independent on each 
element. The formal assembly of the global stiffness matrix and loads is ac- 
complished by introduction of assembly matrices [A^^] through whose use the ele- 
ment level velocity parameters may be expressed as functions of the global 
velocity parameters. For {q} and {6q} we write 
= [AjjHQfQ}: {6qjj} = [A^USQ); 

and from (1.48) thus obtain 

{6Q}'^[Kg]{Q} = {6Q}^{Pg} - {6Q}^[K^]{Q} (1.51) 


In equation (1.51) the global stiffness matrix [K ] and the loads {P_} are 

G G 

defined by 


nelm 

'‘'o' ' 

N=1 

and {Pg} = [Aj^]^{{Fjj> - [0 G]{h“^Pj^}} (1.53) 

The load matrix {P„} contains contributions from the prescribed body for rate 

G 

• 

the relaxation I, and the traction boundary condition The global stiff- 
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ness matrix, as defined by (1.52), will be singular for rigid translations 
(but not for rigid spin). In order to solve the equation (1.51) we define 
a modified global stiffness matrix [K*] and a modified load {P*} as follows: 



p* 

I 


,6lj If (Qj-Qj) or (Qj-Qj) 

, K otherwise 
.Qj if (Qj=Qj) 

NELM _ 

P - Z K^tQt otherwise 


(I. 54) 


(1.55) 


Then (1.51) may be replaced by 

[K*]{Q} = {P*} (1.56) 

If [K*] is not singular, then we solve (1.56) for {Q}, 

{Q} = [K*]"^{P*} (1.57) 

By backsubstitution we obtain the velocity (on the boundary of each element) , 
the spin and the stress rate on each element: 

{qN} = [\][K*]"^{P*} (1.58) 

a 

{g^} = [A^][K*]"^{P*} + {H"^P^} (1.59) 

N 

v(x) = [N(x)][A^][K*]"^{P*> (1.60) 


m(3c^) j 

t(x) - 


QW(2c) 0 
0 QT(x) 


j[H“^G^][A^][K*]"^{p*} + (1.61) 


Equations (1.60) and (1.61) comprise the approximate solution of the boundary 
value probelm. 

Integration of the Motion of the Body ; 

The finite element algorithm just described produces an approximation for 
the stress rate and velocity, as opposed to stress increments and displacement 
Increments. Thus, considerably more freedom of choice of time integration 
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schemes is afforded by the present approach than by Incremental approaches 
(which are predisposed to integration by the relatively inefficient Euler's 
method). In this section we (i) formally state the initial value problem, (ii) 
discuss numerical integration of that problem, and (iii) present a 'forward 
gradient scheme' which stabilizes integration of deformations of bodies which 
exhibit stress relaxation. 

Let {jc} = 2i^»***» be the vector of nodal positions, and let 

{v} = {v^, be the vector of nodal velocities , where ND is the 

12 G • 

total number of nodes. Similarly, let {t} = {t , x ,...,x } and {t} = 

{t^, t^,...,t^} be the quadrature point stresses and stress rates, respectively, 
where G is the total number of quadrature points in the body. To indicate 
the dependence of {v} and {t} on {jf}, {x}, and the time dependent prescribed 
loads, we write* * 

{v} = f[{x}> {t}, t} (1.62) 

{ij = g[{x}, {t}, t} (1.63) 

Since each element node is associated with the same material point 
throughout a deformation, and likewise for each quadrature point, we may write 


each component of {x}, (x), as 

(1.64) 

x^ = (l/jJ)F^ . t^(X^,t) (1.65) 


f = i^(X^,t) (1.66) 

= (l/jJ)F^ . t^(X^,t) (1.67) 


* The function f and g are introduced specifically as a 'shorthand' for the solution 

of the finite element equations, as given by (1.60) and (1.61). In practice 

integrations may be performed on one element at a time. 
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( 1 . 68 ) 


Introduction of (1.64) - (1.67) to (1.62) and (1.63) gives 

(t^> = t] (1.69) 

the definitions of f^ and being clear. Equations (1.68) and (1.69) and 
appropriate Initial values comprise an initial value problem. 

The initial value problem posed by (1.68) and (1.69) and appropriate in- 
itial data is dependent upon the finite element equations. From that same 
discussion, and from the presentation of the finite element equations, it is 
also clear that the finite element-initial value problem is predisposed to 
numerical integration. In this section we indicate the types of numerical 
integration schemes suitable for the present problem, and mention a few im- 
portant differences between the various types. 

The finite element-initial value problem may be Integrated by single step 
explicit schemes, multistep explicit schemes, or (generally multistep) pre- 
dictor-corrector schemes. Three important facts to be kept in mind when 
choosing a particular scheme are 

(1) the solution vector (t„)}, {t (t„) }> at the time t=t„ is of 

N ~T N N 

scalar dimension ND0F+9.G, where NDOF is the number of kinematic 
degrees of freedom of the mesh and G is the total number of quadra- 
ture points. Storage required for implementation or different integ- 
ration schemes can vary appreciably. 

(2) evaluation of <f^, g^> is expensive on account of the complexity 
of the finite element equations. 

(3) the functions f^ and g^ are generally discontinuous at points 
<^, t^> which correspond to material yield surfaces. 

The multistep methods (explicit and predictor-corrector) require rela- 
tively few evaluations of g^^ per step; this is an attractive feature. 
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However, multistep methods are not self starting, the time step is not easily 
changed, they have relatively large storage requirements (since several past 
values of must be carried along), and moreover, they cannot be ex- 

pected to be accurate when the solution crosses a yield surface (since they 
are based on smooth ^polynomial interpolation of the solution over several 
time steps) . 

On the other hand, the single step methods (explicit and predictor- 
corrector) are easily started, the time step size is easily adjusted, and they 
have relatively small storage requirements. They can be expected to perform 
more favorably than the multistep methods when the solution crosses a yield 
surface since smoothing over several time steps is not built in. The disadvantage 
of the single step methods is that a relatively larger number of evaluations of g^^ 

are required per step to achieve a given accuracy when a yield surface is not crossed. 
However, the advantage of single step methods seem to far outweigh the disadvantages. 

In the example accompanying this report the Euler and classical second 
order Runge-Kutta methods were used. Details of these methods may be found 
in many textbooks. Errors of the Euler method were gauged (qualitatively) 
by step-halving and by comparison to results of second order integration for 
randomly selected time steps. Full details are given in the description of 
the example problem. 

It is worthy of special note that complementary work and energy principles 
provide no means whatever for checking the satisfaction of LMB, so it is of 
crucial importance that the numerical Integration scheme not introduce errors 
which tend to unbalance the stress. This maintenance of balanced stress, 
necessary in stress-based finite element algorithms is the counterpart of 
maintenance of compatible deformation, necessary in velocity - based algo- 
rithms. It can be shown that LMB is maintained when the stress t^ is integ- 
rated explicitly, but not when other stresses (such as x) are integrated ex- 


319 



plicitly. Thus, we integrate t^ explicitly, and find x (afterwards) by the 
formula 

T = 1/j F .t (1.70) 

~ T~T -T 

Stability of Numerical Integration of the Initial Value Problem : 

It is possible that the difference between two supposed numerical solu- 
tions of a given initial value problem is much larger than would be expected 
to arise from discretization error alone. As an example, consider integration 
of the stress in a material of the type (1.16) by the Euler method. We suppose, 

3 

for the sake of simplicity, that s(t) is given and f (x) = -2y(-^x ') , so that 
the difference between two solutions satisfies 

^2* “ ~ V(x)]:e(t) - (Syy)^!' (1.71) 

If the elastic matrix and stretching are such that, in the Euclidean norm, 

I I [V(x+Ax) - V(x) ] : e(t) I |/ I I Ax 1 I -»■ 0 (1.72) 

as |(Ax|| 0, then for sufficiently small ||Ax||, equation (1.71) may be 

replaced by 

Ao* = -(Syy)^!' (1.73) 

Defining Ac as Ao = J^-|-Ax' : Ax’ , we may reduce (1.73) to a scalar equation in 
the invariant Aa: 

(Aa) = -(3py)Aa (1.74) 

For an initial value Aa(0) (small), the analytic solution of (1.74) is 

Aa(t) = 40(0)6"^^^^^^^ (1.75) 

Euler's method yields 

40-^ = Aa(0)(l-3yyh)^ (1.76) 

It is clear from (1,75) that Aa decays to zero as time passes. This means 


that the analytic solution of the equation 
6* »= V:e(t) + 1; x(0) * x 

~ z ~ ~ ~o 
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(1.77) 



is stable with respect to sufficiently small perturbations of the deviatoric 
part of T. On the other hand, the numerical solution (1.76) attenuates as 
time passes only if 

1(1 - 3PYh) I < 1 (1.78) 

This means that the numerical solution of (1.77) is stable with respect to small 
perturbations of the deviatoric part of t* only so long as the time step h is 
bounded above as 

|h| < 2/(3 py) (1.79) 

This bound is identical to the bound given by Cormeau [7] (see equations 
16 and 54 in this reference) . 

Time steps such as (W9) are found to be necessary for stability of 
numerical solutions of the finite element -initial value problem presented in 
this report. Argyris et al [8] remark that this time step restriction amounts 
to limiting the inelastic strain increment to be smaller than the elastic strain. 
Since the elastic strain is usually very small in metals such as those used in 
structures, this implies that a finite deformation analysis would entail an 
impractically large number of steps. 

The work of Kanchi et al [9] and Atluri and Murakawa [4] suggests the 
modification we now describe. To Improve the estimate of the Inelastic strain 
Increment in a time step, we replace by estimate of the mean value 

of the Inelastic stretching in that time step: 


de^^ 

eP(T(tjj+0h)) = e*’(T(t^)) + 6h 


: 6 * 


I=In 


(1.80) 


where the parameter 6, O_5.0^1, serves to locate the time at which the mean 
value is achieved. As 6 goes from zero to one, the estimate of the creep- 
stretch becomes increasingly more conservative. 
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Equation (1.80) may be introduced to the finite element algorithm 


through the constitutive equation; (1.16) becomes 


“ Ye-5 ?e 


where 


Yj - [f ^ 1 ; Eg = -VgicP 


de 


P -1 


From V. (1.81) we derive W. just as we derived W from V: 

:: y z z 

te " Ye ■ I’ '^ijkl " 2 ^’^ik'^lj '^ik'^lj^’ 


= -(W^+T):eP 


=0 ~e ~ 


‘~e 


(1.81) 


(1.82) 


When a material which exhibits relaxation is to be analyzed, Wg and Eg are in- 
troduced to the finite element algorithm for W and E. 

Example ; Growth of a Void in a Viscoplastic Medium ; 

In this example we examine the growth of a void in a hypoelastic/visco- 
plastlc medium. This problem has been studied (numerically) by Burke and Nix 
[10], who treated the material as rigid/viscoplastic. We present the problem 
as a demonstration of the performance of the finite element algorithim. The 
material exhibits stress relaxation, so the forward gradient scheme must be 
used to stabilize the time integration. The present results agree quite closely 
with those of Burke and Nix. 

The motion is assumed to be plane strain, and throughout the body is a 
doubly periodic array of cylinderical voids. Due to the symmetry we need analyse 
only one quadrant of one rectangular cell of the body. The finite element 
mesh and boundary conditions are described in Figure 1. 

Burke and Nix motivate their study by explaining that certain theories 
for the initiation of creep fracture supose that the growth of voids can be 
"attributed to the inhomogeneous plastic deformation of the surrouding grains." 
Furthermore, "finite fracture strains can be predicted only when a void lies 
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In the neighborhood of another void." Such a study necessarily involves a 
number of special cases. For our purposes, that of demonstration, only one 
case is taken. 

The problem has been analysed in three parts. In the first part the 
cell is brought rapidly from the virgin state (stress-free) to a state of 
purely elastic strain. This is accomplished by a single RK2 step. In the 
second part, relatively small time steps are taken while the stress relaxes 
from the elastic distribution to a nearly steady creep distribution. In the 
third part, time steps are taken which produce 1% nominal elongation of the 
cell in each step. To stabilize time integration in the second and third 
parts the forward gradient scheme is used, the stability parameter 0 set as 
0 = 1/2 and 3/4. Only the Euler time stepping scheme has been used in the 
second and third parts of the problem. 

The material model is a special case of (1.16); 

® j. P 
e = e + e 


= (^)g* - (^)(I;g*)I 



This model corresponds to that of Burke and Nix with (their) creep exponent 

-19 -1 

n=l. The fluidity y is set as (psi-sec) . The velocity at the 

top of the cell (see Figure 1) was adjusted so that a specimen with no void 

would experience a homogeneous constant stretching e^^of e^ =0.25x10 sec.”^. 

Since the material was treated as rigid/viscoplastic by Burke and Nix, our 
choice of elastic constants is somewhat arbitrary. We have taken Young's 
modulus E=3xl0^ psi and Poisson ratio V"0.4, so the material is somewhat like 
steel in its elastic response. 
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11 33 

In Figures 2, 3, and 4 the contours of stress t , stress t and mean 

stress, have been plotted for L (the elongation of the cell) L=1.01. The 

3 

stress concentration where the hole edge crosses the x axis is approximately 
2.7*. This is quite reasonable since the theoretical value for an isolated 
void in a purely elastic medium is 3.0. In Reference 10 an approximate value 
of 2.66 was found for the rigid plastic material. In Figure 5 the contours of 
effective strain are plotted for L=1.01. Qualitatively this com- 

pares very well to Figure 7 in [10]. 

In Figure 6 the deformation is traced from L=1.0 to L=1.5. These de- 
formations are physically tenable. We remark that no indication of any num- 
erical instability was observed in the course of integrating thisi deformation. 

11 33 

In Figures 7, 8, and 9 the contours of stress f ^ l ^ and mean stress, 
have been plotted for L=1.50. They compare very well to the stresses found 
in [10] (see Figure 8 there). We note that the stress concentration has drop- 
ped to 1.71. The stress concentration depends strongly on the geometry of the 
specimen; as such, it was observed to decline steadily throughout the deforma- 
tion. In Figure 10 the contours of effective strain rate are plotted for L=1.5. 
Again, the qualitative agreement with the results of Burke and Nix [10] is noted 
(see Figure 9 there) . 

We conclude by noting that in the present analysis only 56 four noded 
elements were used, as compared to 56 eight noded elements used in the analysis 
of Burke and Nix. Considering the agreement between their results and our own, 
the present method appears to have performed very well, in spite of the large 
disparity in the degrees of freedom of the finite element mesh. 


* A stress concentration of approximately 2.59 was observed for the elastically 
stressed medium. 
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Part II ; Fracture Analyses Under Creep 

Numerous experimental studies have been undertaken with the purpose of 
finding a parameter which correlates with creep crack propagation rate. Most 
of these investigations consider as candidate parameters, K^, some form of net 
section (or reference) stress or in more recent studies C*. See, for example, 
the review article [11] and [12-14]. Since the introduction of C*, there appears 
to be less emphasis on as a parameter, however, there are apparently real 
materials and conditions for which either net section stress or provide better 
correlation with crack growth rate than C*. 

As illustrated in Fig. 11, the above three parameters might be expected 
to correlate three distinctly different creep crack growth situations. In Fig. 
11a, a crack and its associated ligament are shown for a material and geometry 
which results in negligible creep strains everywhere except in the vicinity 
of the crack-tip. This condition is analogous to that of small scale yielding 
in elastic-plastic fracture. 

Fig. 11b represents a situation in which. C* might be considered an appropri- 
ate parameter. This situation is characterized (i) by the body being essentially 
at steady-state creep conditions (which implies very slow crack propagation) 
and (ii) by the creep-damage process-zone being local to and therefore controlled 
by the crack-tip field. Fig. 11c illustrates the type of situation for which 
net section stress might be expected to control crack growth. In this case, 
the main feature is the widespread creep damage zone. 

It is seen from Fig. 11 that Intermediate situations can occur. For 
example, suppose a particular material and geometry results in a crack propa- 
gation rate such that elastic strain rates are not negligible compared to 
creep rates (i.e., non-steady creep) and at the same time, creep strains are 
no longer localized to the crack-tip region. While neither nor C* could 
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be valid parameters for this case, it appears reasonable to expect that crack 
growth rate is still determined by the local crack-tip field since the creep 
damage process zone is still assumed to be local to the crack-tip. 

In the present study, we are concerned primarily with behavior bounded 
by at Illustrated in Fig. 11a and Fig. 11b. That is, we consider conditions 
in which the creep damage zone and presumably crack propagation speed are con- 
trolled by the crack-tip field. Therefore, if we have a parameter which 
characterizes the crack-tip fields during such behavior we presumably have a 

parameter which will characterize creep crack propagation rate. A parameter 
which spans the gap between ICj. controlled growth and C* controlled growth 
has been introduced in [15] and subjected to initial scrutiny in [16]. This 
parameter is referred to as defined by a path-independent, vector 

integral. For stationary cracks, it has been shown [15,16] that the related 
quantity a measure of the amplitude of the HER crack-tip field which 

presumably exists for both non-steady and steady-state creep. It has also 

• • dU 

been shown that has the energy Interpretation ~ ^ non-steady 

as well as steady-state creep. 

In the process of exploring this new parameter, it has been found [16] 
that despite C* being a valid crack-tip parameter for strictly steady-state 
creep conditions, it is not equivalent to the parameter under any 

conditions and therefore does not have the energy interpretation commonly at- 
tributed to it. Since experimentalists use the energy interpretation as a 
means of "measuring" C* it seems more appropriate to refer to these experi- 
mental results as (T, ) . 

1 c 

To start with, we define ^ and a generalized vector integral 
and also summarize the properties and the relationship of these parameters. 
The remainder of the paper discusses several finite element calculations for 
both stationary cracks and propagating cracks. The crack propagation study 
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uses a combination of analytical, numerical and experimental results to show 


that creep crack growth in 304 stainless steel at 650 C occurs under essen- 
tially steady-state creep conditions. Finally, based on this observation, a 
simple crack growth prediction methodology is outline. 

Constitutive Equations : 

In this study we assume strains are infinitesimal and the deformations 
small. Furthermore, we assume material behavior of the type: 


- .e,.c _ . |3. X n 1 I 

ij ij ij ijk£ 2 ' ij 


(II. 1) 


• 0 • c 

where and are the elastic and creep strain rates, respectively, 

is the tensor of elastic moduli, is the stress rate, is the deviatoric 

stress (Tl,==t.,- and a is the equivalent stress given by 

ij ij 3 kk ij 

- 1/2 

a = (3/2) (t^jT . The parametersy and n are those of the familiar Norton 

power law: 


€ = Y (a)^ 


(II. 2) 


where 


€ = [(2/3)c..I.^]^/2 

The constitutive law (H-D can result in steady-state creep response (i.e., 

after some period of time provided the boundary conditions are some combination 

of time invariant tractions or time invariant displacement rates. 

Fracture Parameters (AT) and C*: 

— c — 

We now define two vector quantities which have applications to fracture 
analysis under creep conditions. The first quantity is (&X)(, recently de- 

fined by Atluri [15] and subsequently examined in greater detail in [16]. In 

[15] (at) is defined in the context of finite strains and large deformations. 

— c 

Here we give the corresponding definition for infinitesimal strains and small 
deformations. 
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^ “o 


/ [n,AW-n. 

Jj- i J 


3Au, 


- / [n^AW-n,(T.^+it.^) — 

^ 234 ^ 


•]dS 


Lt 

£•>0 


X-V " 3x^ 


y* n^AWdS + J n^AWdS - ^ 
*^^12 *^^45 


9A 


\ 


'k 3 X, 


dS 


-X 


45 

3Au, 


n (t .,+At,, ) -T dS 

j jk jk 9x. 


(II. 3) 


The various contour integral paths and their outward unit normals n as well 

as V and are illustrated in Fig. 12 for a two-dimensional, cracked body. 

In writing (II. 3) it has been assumed that S +S^=r, -+r,c- where S and are the 
^ e t 12^ 45 e t 

portions of the crack surfaces with applied Incremental displacements, 4Uj^, 
and applied tractions t, , respectively- The initial stress for the Increment 

tC 

is denoted The mass density is P and the acceleration and body force 

components at the end of the increment are a^ and fj^, respectively. The 
quantity AW is the Incremental stress-working density and is given by 


AW = (t . . + -r At . . )Ae . . 

13 2 xj X3 


(II. 4) 


The right equality of (II. 3) shows that (4T^)^ is independent of the selection 
of 1234 the fields within V-V^ are sufficiently well behaved for the 
divergence theorem to be applicable) . It is important to note that this path- 
independence exists during non-steady as well as steady-state creep. 

In the present study we consider cracks along tte x^ axis and symmetrical. 
Mode I type deformations. Furthermore, we consider traction-free crack sur- 
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faces and assume body forces and accelerations are negligible. Under these 


conditions, only is of interest and we have 


- r 

(AT ) = / [n AW - n.(T.,+AT.,) -r ]dS - / - ^ Ae .. dV (II. 5) 

1 c / ^1 3 jh jk 3x J 3x 3 k 

‘234 X V X 


where we have now taken the limit of the volume integral. 


It has been shown in [15] that kas the physical meaning* 

.AU -AU 


(II. 6) 


where AU2 and AU^ are the incremental potential energies for two cracked 
bodies which are indentical in loading history and geometry except that the 
second body has an incrementally longer crack by the amount dc^^. In creep 
applications, it is convenient to define the quantity 


(Tj 


_ Lt ^^"^l^c - ^^^l^c 


(II. 7) 


1 c At-»-0 At At 

where At is the time increment. Comparing (II. 7) and (II. 6) it can be seen that 

(T ) has the physical meaning which is commonly attributed to C*, i.e. 

X C X 


We now state a generalized definition for the C* parameter which has 
been derived in [16]. 


*The sign convention for AU and AU is reversed from that of [15,16] to reflect 
the conventional definition of potential energy. 
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(II. 9) 


Lt r ^ \ 

= In / [n.W*-n.T ^]dS 

1 £->0 1 j jk 3x^ 

^ e 

■ £ K«*-Vjk 

*^^234 ^ 

+ ^ n I r P(a,“fi,) ^ dV + r n W*dS + /* n i 

^^0 (Jv-v/ ^ ^ Jr^2 *'^45 

^ 9 u ^ 9 u I 

- / t, ^ dS - / n.T., ^dsl 

Js k 9x. Jg j Jk 9x. J 

t e ' 


W*dS 


Based on the same simplifying conditions used in obtaining (II. 5) we have 


r 

C* = / [n W*-n,T T— 

1 Jr„, 1 J J’' 


(II. 10) 


where it is seen that the volume integral no longer is present. 

The quantity W* which appears in (II. 9) and (II. 10) is usually defined as 


t . . 

riJ 

= / T..de.. 

J Q X3 X3 


(II. 11) 


Using the steady-state case of (II. 1) and the associated incompressibility condi- 
tion, the following more useful expressions can be derived [ifi]; 


w* • ^ “ 


w* - 


(11.12) 


(11.13) 


As noted previously, C* is often stated to have the energy interpretation 

which was given for (T.) in (II. 8). It has been shown in [16] that this Incor- 

1 c 

rect. The relationship of the steady-state value of (T.) (i.e.(T, ) ) and 

■ 1 c 1 css 


C* is given in [16] as: 


^^l^css *^1 ^n+1^ e-»-0 


J" n^(a)'^^^dS 


(11.14) 
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Approximate numerical evaluation of (H.l^) in [16] has shown that 

C* agree to within 2% for plane strain and differ by as much as 14% for plane 

stress. 

From the above discussion it is clear that C* and (T) are not equivalent 

— — c 

quantities under any condition despite their being derivable from the same 
conservation law. The quantity (X)^. follows more directly from the conservation 
law and is the more general quantity not only in that it is applicable to non- 
steady as well as steady-state creep but also in that it is applicable to 
constitutive laws which are more general than (II. 1) [15]. The quantity relies 
on the special property of (II. 1) which allows the existance of a potential W* 

for the deviatoric stresses, Furthermore, since W* does not have a 

« 

physical meaning whereas W has the meaning of rate of stress-working density, 

it is understandable that energy interpretation whereas does 

not. In light of this conclusion, it seems more appropriate to refer to experi- 

mental measurements of " as measurements of (T-) as opposed to measure- 

X C 

ments of C* or etc. 

Finite Element Equations ; 

The following summarizes the finite element model. For a more complete 
description see [16]. The model is based on the principle of virtual work: 


f x..6e..dV - / t.Su.dS = 0 

Jy 13 13 L 1 1 


(11.15) 


By substituting the following incremental stress-strain relation 


{t>^ = + tE]{Ae - [E]{Ae^}. 


(11.16) 


into (15) and applying customary procedures we have the final equation: 

[K]{AQ>^ = {T}_ + {S - {R},. . * 

i i cl 1—1 


(11.17) 
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{R}j_j^=2: r j_^dV (11.21) 

e Jv 

e 

The form of (11*17) makes this an Initial strain formulation. It can be 

seen that [K] is the elastic stiffness matrix and remains xinchanged through- 
out the time Incrementing process. The quantities {Ae^}^ in (20) are predicted 
prior to the solution of (H*17) using (II. 1) and Having sovled (11.17), and 

thus obtained {Ae}^, the actual values of and {AtI^ are obtained by 

subdividing the time step and performing an Eulerian Integration based on sub- 
increments of {Ae)j. As a result of this Integration procedure, better ad- 
herence to the postulated constitutive law (II. 1) is achieved but at the price of 
introducing some disequilibrium (i.e. {R}j?*{T}j). This disequilibrium is cor- 
rected, however, in the next time step as a consequence of {R}^. ^ appearing 
In (11.17). 

The time step size for the calculations is automatically regulated based 
on two criteria. The first criteria is the maximum error in the predicted 
creep strain Increments used in solving (11.17) as compared to the creep strain 
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increments from the subsequent integration procedure. The second criteria 
is the maximum creep strain increment compared to the total elastic strain. 

In the present study, the criterion for maximum error in predicted creep 
strain is 20% and the criterion for maximum creep strain increment is 100%. 

For problems which have been considered, it appears that the above model and 
criteria give accurate transient solutions and converged steady-state solutions 
with time step size comparable to those used with more expensive tangent 
stiffness methods. 

Verification of Model ; 

A compact specimen has been chosen for verifying the model. The parti- 
cular geometry and materials were chosen to coincide with those used by 
Ehlers and Riedel [17] and are Illustrated in Fig. 13 along with the two finite 

element meshes used in the verification. Both meshes consist entirely of eight- 

noded isoparametric elements, assume plane strain conditions and use collapsed 

(i.e. triangular) elements at the crack-tip. For the 102 element model these 

- 1/2 

crack-tip elements are given a singular strain field (r ) by shifting the 
appropriate midside nodes to their quarter-points. The 300 element model uses 
a non-singular crack- tip. 

The elastic for the 300 and 102 element meshes are 24.1 and 24.3 N/mm, 
respectively and agree with the value 24.2 N/mm from Srawley [18] to well within 
1%. The steady-state (t=600hr) values of C* for the 300 and 102 element meshes 
are 131 and 130 N/m*hr, respectively, and" agree well with 134 N/m*hr from Shih 
and Kumar [19] and 137 N/m.hr (t 300hr) from Ehlers and Riedel [17]. Based 

on these results it is concluded that the numerical procedure in general and 
the quarter-point crack-tip elements in particular, are accurate and efficient 
tools for creep fracture analysis. 
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Calculations of (T,) and (T, ) : 

1 c 1 c 

Now we consider calculations of (T. ) (i. e. (AT,) /At). The 300 element 

1 c 1 c 

mesh results for (T ) , as computed from (II. 5), are shown in Fig. 14 as the solid 

1 c 

« 

curve. This curve shows the time dependence of (T^) during the non-steady 

c 

portion of the creep calculation. The steady-state value of is 130 N/m'hr 

and thus is in agreement with the previously mentioned relationship between 

(T ) and C*. In [16] it was found that the evaluation of (II. 5) using the 102 
1 css 1 

element mesh gave values of (1^)^, which were generally in poor agreement with 
those of the 300 element mesh. The volume integral of (H*5) was determined to 

be the cause of this behavior and it was supposed that the origin of the 

- 1/2 

problem was the use of the r strain singularity as opposed to the HRR type 
singularity (i.e. r However, several calculations with special 

conforming elements which Impose the HRR type radial dependence of strain [20], 

have shown that this is not the case. 

It now seems that the difficulty experienced in computing (T ) when 

1 c 

using singular crack- tip elements is related to the existence argument for the 
limit of the volume integral in (II. 3). For the case when the asymptotic field 
has singular radial dependence but does not identically satisfy the following 
condition on angular behavior. 


Lt y 9T^i-^( g»Q) 
£+0 


L 


9x, 


Ae..(e,0)d0 = 0 


( 11 . 22 ) 


the subject limit in (11,3) does not exist (as discussed in Appendix A of [16]). 
The condition (11.22) need not be satisfied exactly if one does not ahve a singu- 
lar radial dependence as is demonstrated by the results from the 300 element 
mesh. 

The efficiency, simplicity and general accuracy of the quarter-point ele- 
ment procedure makes it a very attractive alternative to the use of very re- 
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fined nonsingular meshes or the derivation of singular crack-tip elements which 
satisfy (11.22) a priori. Therefore, a practical solution to this problem is 
sought. As noted above, the difficulty is associated with the volume integral 
over the singular elements. Therefore, calculations were made in which the 

volume Integral over the crack-tip elements was omitted. The resulting quantity, 

♦ 6 

which we call (T-) , can be written 
i c 





[n. AW-n . (t +At ) 
1 J Jk jk 




Ae,, dV 


(11.23) 


where consists of the singular crack-tip elements. The dashed curve of 
Fig. lA is from the 102 element mesh. It can be seen that coin- 

cides with the solid curve for times after about 30 hours. For this mesh 

and problem it can therefore be said that (T,)^ is a valid path-independent, 

1 c 

crack-tip parameter for times after 30 hours and for values of beginning 

at approximately 1.6 of the steady-state value. The steady-state parameter 
is still significantly path-dependent at 30 hours. 

For the 102 element mesh, the crack- tip elements are 5% of the ligament 

• 6 

size. We therefore assign 5 a value of 0.05. A quantity similar to (^ 2 ^)^, 
was computed using the 300 element mesh. In this case, a semi-circular 
region of radius approximately 3% of the ligament was omitted from the evalua- 
tion of the volume integral of (H*5). This result which we deonte 

e=3% is also shown in Fig. 14. This curve seems to indicate that the validity 

* 5 

of (T|^)^ can be expanded to earlier times with rather moderate reductions in 

the crack-tip, quarter-point element size. For example, the results of Fig. 14 

• 6 

indicate that a 5 of 3% of the ligament would result in being valid as 

early as seven hours and for as large as 4.3 its steady-state value. 


Creep Crack Growth in a Strip : 

We now consider the problem of a finite height (2h) infinitely wide strip, 
with a serai- infinite crack. Loading consists of uniformly applied displace- 
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ment rates (6) at the top and bottom edges (y=+h) such that Mode I behavior 
results. This problem has been chosen for two reasons. First, since the strip 
is infinitely wide and the boundary conditions do not change with time, the 
propagating crack-tip fields can be expected to reach a "convecting steady- 
state" creep condition. Here we use the phrase "convecting steady-state" 
to mean that the field remains unchanged in time with respect to a coordinate 
system which is centered at and moving with the crack-tip. This terminology 
is used so as not to confuse this condition with the usual steady-state creep 
condition in which material stress rates are zero. 

The second reason for choosing this problem is that C* can be evaluated 

analytically for the special case of steady-state creep (stationary crack). 

The analytical evaluation of (II. LO) follows easily if one chooses a rectangular 

contour in which the horizontal portions coincide with the top and bottom 

edges of the strip (i.e. y=+h) and the vertical portions are at x=^«. In such 

a contour one finds only the vertical portion at xC=+» is non-zero and therefore 

C* = 2hW* (11-24) 

1 «> 

For the corresponding elastic problem with applied displacement 5 , one finds 
a similar relation. 

J. = 2hW (11.25) 

X 00 

It has been noted that (T, ) and C* are related and therefore it is possible 

1 css 1 ^ ^ 

to obtain (tj from (H-14) and (11.24). The direct evaluation of (T^^) in terns 
1 css 

of either its integral representation (H-5) or its energy representation II. 6) 
requires knowledge of the stresses in the region of the strip adjacent to the 
crack-tip and therefore is not a trivial task. 

The material properties used in this problem are representative of 304 
stainless steel at 650°C. These material properties and the finite element 
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Table 1. Summary of Analysis Parameters for Creep Crack Growth 
in the Plane Strain Strip of Fig. 15 


from (6.4) 

Analytical Results Computed Results (mm/hr) 


c* 

1 

(N/mm'hr) 

remote 

T 

yy 

(MPa) 

6 

(mm/hr) 

5 

(elastic) 

(mm) 

•^1 

(N/mm) 

(N/mm* hr) 

•^1 

(N/mm) 

average 

upper 

bound 

0.05 

83 

3.44 X 10"^ 

5.04 X lO"^ 

4.18 

4.99 X 10~^ 

4.19 

1.00 X lO"^ 

5.00 X lO"^ 

5.0 

148 

1.94 X 10"^ 

8.95 X 10“^ 

13.2 

4.99 

13.2 

2.22 X lO"^ 

1.11 X lO"^ 

50 

197 

1.45 X 10“^ 

1.19 X lO"^ 

23.5 

49.8 

23.5 

3.30 X 10“^ 

1.65 


discretization are given in Fig. (15) • Note that collapsed, eight— noded, quarter- 
point elements are used at the crack-tip. 

The mesh for this problem may at first appear rather coarse; however, 
elastic and steady-state creep solutions obtained with this mesh are suffic- 
iently accurate to justify its use for the study at hand. The comparison of 
computed elastic values and steady-state C* values with their analytic values 
is given in Table 1. 

The first step in this numerical study is to select three values of C* 
which span the range of values reported in the literature for 304 stainless 
steel at 650°C.f The values which have been chosen are 0.05, 5.0 and 50.0 N/mm'hr. 
Having these values, the remote (x=“) steady-state are determined as well 
as the edge displacement which results in the same remote elastic '^yy* 

These displacements are applied to the model elastically at t*0. Next, the 
steady-state edge displacement rates are determined analytically. Using the 
elastic solution as an initial state, the displacement rate, 6, is applied 
until the model reaches steady-state. 

The next step in this study involves the selection of upper bound crack 
velocities for the three chosen values of C*. The following formula is based 
on the experimental data reported in [13,14] and represents data from center- 
crack, double-edge-crack, compact, and round bar specimen types. 




where a 


1.68 X 10 
3.36 X 10 


(upper bound) 
(average) 


(11.26) 


tThe use of C* rather than (T.) is due to the existence of the analytical ex- 
1 i c 

pression (24) and is justified by the numerical similarity to plane 

strain conditions. 
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Having reached steady-state, the crack is propagated at the upper bound 
velocity of (11.26) until it is determined that a convecting steady-state has 

been reached. 

The crack growth simulation is accomplished through a combination of 
mesh shifting and periodic remeshing as illustrated in Fig. (16) . The region A 
represents the quarter-point elements which remain centered about the crack- 
tip. The B type elements are standard eight-noded isoparametric elements which 
distort during mesh shifting so as to keep the region A centered at the crack- 
tip. The procedure is to shift the region A (and thus the crack-tip) by 
shifting appropriate nodes of the region A and type B elements. This shifting 
is done without altering element connectivity. Eventually the type B elements 
become overly distorted at which time the element connectivities are redefined 
in the vicinity of the crack-tip so that additional shifting is possible. 

Each occurance of shifting or remeshing requires that shifted nodes have 
their displacements interpolated and that shifted elements have their 2x2 
Gauss point stresses interpolated. The displacement interpolation is by the 
usual isoparametric shape functions. The stress interpolation uses linear, 
two-dimensional Lagrangian polynomials in element local coordinates. In the 
following calculations, the nominal size of the crack growth increments is 0.4 
mm or 2% of the crack-tip element width. For the highest velocity case 
(C*=50 N/mm*hr), this results in crack growth at approximately every fifth 
solution step. 

Results for a Plane Strain Strip ; 

The results of the plane strain strip calculation with C*=50 N/mm*hr 
and =1.65 mm/hr are given in Fig. 17. The values of (Tj^)^ C* are given 
for the portion of the calculation prior to steady-state as well as during the 
crack propagation portion. The band represents the range of values obtained 
from the four contours Illustrated in Fig. 15. Both converge 
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to the 50 N/mm'hr value at steady-state. During the crack propagation, it is 

seen that and C* do not depart significantly from their steady-state 

value. This means that this combination of loading and crack -speed results 

in the crack-tip fields being essentially at steady-state conditions. This 

* 5 * 

in turn means that both (or and C* are valid crack-tip field 

parameters during crack growth. 

A closer view of the crack propagation portion of these curves is given 
in Fig. 18. The dashed curves bracketing the initial portion of the solid 

curves represent the degree of path-independence and continue to be represen- 
tative of the path-independence observed during the crack propagation steps. 

• 6 

For both (1^)^ and C*, it is seen that the strip has essentially returned 

to its steady-state condition prior to each crack growth Increment. It is 

• 5 

thought that the larger departure of (l^^c from steady-state (as compared to 
C*) is more representative of the non-steadiness of the crack-tip field since 
the validity of C* in general and the numerical evaluation of W* (11.13) in par- 
ticular, are based on the existence of steady-state conditions. 

The effect of remeshing is seen at approximately eight hours. The first 


two steps after the remeshing were found to result in rather erratic contour 
integral values and are not indicated in these figures. The equilibrium cor- 
rection feature of the present model and the automatic time step regulation 


procedure both act to quickly restore equilibrium at the crack-tip. 

The propagation portion of the calculation with C*=5 N/mm*hr and =0.111 
mm/hr is given in Fig. 19. Here again it is seen that both (T^)^ and C* have 

converged to the analytical value of C* (to within two percent, which is also 
about the degree of path-independence) . Comparing these results with those in 
Fig. 18 for the higher C* and crack speed it is seen that steady-state creep 

conditions were not reached until 12 hours as opposed to approximately two 


340 



hours in previous cases. Also, the return to the steady-state value after 
mesh shifting takes more time (two hours compared to 0.25 hours). However, 
when compared to the time between crack growth steps (both use-0.4 mm) it is 
seen that the lower velocity case returns to steady-state well before the next 
growth step occurs. This result indicates that lower load levels and crack 
speeds are inherently closer to steady-state conditions. While this behavior 

seems intuitively correct, it should be kept in mind that these results depend 
on the empirical formula (H.26) which is valid only for 304 stainless steel. It 

remains to be seen if similar behavior occurs in other materials. 

A calculation has also been done for the case of C*=0.05 N/ram*hr. As 
a result of the large number of solution steps between crack growth steps, when 
using the maximum velocity of 5 x 10 ^ mm/hr, the calculations used a higher 

_3 

velocity (5 x 10 mm/hr). Even at this \mrealistically high velocity (for 
this level of loading) , the behavior is more steady-state-like than the 
case of C*=5.0 N/mm.hr described above. 

Creep Crack Growth in Double-Edge-Crack Specimens ; 

The purpose of considering this problem is to apply the model to a 
problem for which experimental data exists. While much experimental data has 
been reported in the literature, most authors do not include sufficient in- 
formation to allow a numerical simulation of their experiments. The current 
problem is based on the experiments of Koterazawa and Iwatwa [21]. The primary 
reasons for selecting this work for study are the crack length versus time 
histories were given and that the experiments were performed on 304 stainless 
steel for which high temperature elastic and creep properties were already 
available. 

The geometry of the experimental specimens is given in Fig; 20. The 
finite element mesh for the calculations is shown in Fig. 21 with contour 
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integral paths being indicated by dashed lines. It can be seen that the mesh 
takes advantage of the two planes of symmetry for the specimen and does not 
model the 60° notch. The initial crack length indicated in Fig. 21 corresponds 
to the notch depth in the specimen. All calculations for this specimen assume 
plane stress conditions and use the material properties given in Fig. 15. 

Elastic results for two crack lengths are compared in Table 2 with those 
(based on formulas for K^) from [22] and are seen to be in good agreement. 

The material properties are those of 304 stainless steel at 650°C and 
are assumed to be the same as those used in the strip analyses. (See Fig. 

15). Calculations have been made for remote applied stresses of 157 and 176 MPa. 
The experimental crack growth histories for these two stress levels are repro- 
duced from [21] in Fig. 22. It is seen from these curves that the first two- 
thirds of the specimen lives are characterized by crack velocities of less than 
0.01 mm/hr compared to nearly 0.5 mm/hr as rupture is approached. 

The primary purpose of the following calculations is to verify the 
conclusions which were reached in the previously described strip calculations; 
that is, that the crack-tip fields are essentially creep-steady fields 
even for the most rapid creep crack velocities. These calculations will be a 
valid check because the input to the calculations is only the remote applied 
stress and the measured crack velocity history, and does not in any way depend 
on experimental determination of C* or did the strip calculations. In 

fact, Koterazawa and Iwata do not report such measurements in [21]. 

Analysis of Initial, Low Velocity Crack Growth ; 

This section describes the simulation of the initial portion of the 
crack velocity histories given in Fig. 22. In all of these calculations, the 
entire load is applied elastically at t = 0 and held constant throughout the 
subsequent creep solution steps. The convergence of and C* to their 
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steady-state values is shown in Fig. 23, with the dashed lines in the C* 
plots denoting the degree of path-independentce. It is seen that steady-state 
conditions are reached between a half and one hour after the load is applied. 

(Table 2 summarizes the computational aspects of this portion of the calcu- 
lation.) Therefore, it is seen by refering to Fig. 22 that crack growth does 
not begin in the two specimens until well after steady-state conditions are 
reachdd. Since the current calculations assume small displacements and infinitesi- 
mal strains, and since only the strain and displacement magnitudes depend on time 
once steady-state is reached, there is no reason to continue the numerical calcu- 
lations to the crack initiation times indicated by the experimental results. 
Therefore, the initial crack propagation is simulated at times after steady-state 
conditions are reached but much earlier indicated by the experiments. 

The crack growth simulation results are shown in Fig. 23. The crack 
increment size for this study was approximately 0.01 mm which is nominally 2.4 

percent of the crack-tip element size. It can be seen that only one mesh shift 
(l.e., crack growth step) was modeled. It is clear from this figure that the 
time it takes for the specimen to return to steady-state is significantly less 
than the time to the next crack growth increment (indicated by dashed lines). 
Therefore, the initial portion of the crack growth histories of Fig. 22 are 
clearly occurring under essentially steady-state conditions and thus C* as 
well as (T^)^ are valid crack-tip parameters. Since an increase in C* re- 
sults in a more rapid return to steady-state conditions, the above conclusion 
will remain valid for the initial constant velocity portions of the curves of 
Fig. 22. 

When crack growth occurs so slowly that the crack-tip is essentially 
at steady-state, the crack-tip field does not depend on the history of the 
specimen. Therefore, assuming steady-state conditions continue to exist, it is 
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possible to skip to the final stages of crack growth without modeling the 
intermediate crack growth. If it is found that crack growth is still slow 
enough for steady-state conditions to exist, then it seems reasonable to expect 
that the behavior at intermediate crack lengths is also of a steady-state type. 
The following describes the results of this procedure when applied to the 
two double-edge-crack specimens. 

Analysis of Final Stage of Crack Growth ; 

To analyze the final stage of crack growth, the crack length is in- 
creased to 2.75 mm and the process of applying the load elastically and creeping 

to steady-state is repeated. Table 2 summarizes the computational aspects of 

. ,5 

this process. The convergence of and C* is their steady-state values is 

shown in Fig. 24. Having reached steady-state, the cracks are grown at the 

rate suggested by the last portion of the crack histories (Fig. 22) as shown in 

Fig. 24 . The significant increase in the frequency of mesh shifting (compared 

to that in Fig. 23 due to the velocity increase makes the details of the 

curve difficult to distinguish in this figure. However, the time step size is 

such that six or seven steps occur between each crack growth increment. Unlike 

. ,s 

the strip problem, the values of and C* are clearly increasing during 

this crack propagation process. 

It is necessary to determine whether this Increase in the crack-tip 
parameters is due to the crack-tip no longer being at steady-state conditions 
or whether it is due to the crack-tip no longer being at steady-state conditions 
or whether it is due to the incraese.in crack length. This is accomplished by 
continuing the calculation without further crack extension. If the value of 
the parameters do not change significantly with time, this means the increase 
was largely due to the crack length Increase and that crack growth is still 
occurlng under essentially steady-state conditions. Examination of the final 
portions of the curves of Fig. 24 shows that this is the case. 
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Table 2. 


Computational Aspects of the Elastic and Non-Steady Creep 
Portion of the Double-Edge-Crack Calculations 


applied 

stress 

(MPa) 

crack 
length 
( mm) 

Elastic Solution 

Creep Solution 

total 
CP time 
(sec) 

Ji 

( N/mm) 

difference 
from 111 

W 

CP 

time* 
( sec) 

At (hr) 

initial/final 

steps to 
steady-state^ 

157 

1.75 

1.12 

(- 2 . 1 ) 

38 

8 • 10 " Vl.9 • 10 "^ 

90 

795 

176 

1.75 

1 .40 

(- 2 . 6 ) 

38 

8 • 10"®/^. 5 ■ 10 "^ 

100 

880 

157 

2.75 

1.79 

(- 3 . 3 ) 

38 

4 • 10 “®/ 8.6 • 10 “^ 

211 

1820 

176 

2.75 

2.25 

(- 3 . 2 ) 

38 

2 • 10 "®/ 4.4 . 10 "^ 

205 

1770 


^Control Data CYBER ?4 

^Solutions are stopped at times indicated in Figs. 5,19 and 5.20 



Based on this analysis; it appears that the conclusions reached as a 
result of the strip calculations are still valid. Since, (1) the strip analy- 
ses are much less expensive than this analysis of the double-edge-crack geometry, 
(il) the steady-state C* for the strip is easily obtained analytically and (ill) 
the crack-tip parameters do not depend on crack length for the strip geometry, 
it seems that similar studies for other materials and/or other temperatures could 
most effectively be accomplished through the use of the strip geometry. The 
need for such studies follows from the vast simplification of fracture analysis 
and prediction which results if crack growth occurs under steady-state conditions. 
More will be said about this point in the conclusions. 

Summary and Conclusions ; 

It has been noted that despite the fact that C* characterizes the crack- 
tip field under steady-state creep conditions, it does not have an energy or 
energy rate interpretation. A related path- independence integral parameter 

(T^)^, however, does have the energy rate interpretation commonly attributed 
to C*. The derivation of (T.) does not rely on the existence of steady-state 
creep conditions and thus is a valid crack-tip parameter for non-steady creep 
conditions as well as for steady-state creep. 

An Initial strain finite element approach provides for improved 
adherence to postulated constitutive behavior and for equilibrium correction 
has been summarized. The accuracy and efficiency of this model with eight- 
node isoparametric elements and the quarter-point crack-tip element approach 
have been verified through several calculations for a compact specimen geometry 
and a strip geometry. Also, a method of simulating crack growth through 
shifting of the quarter-point singularity elements and periodic remeshing has 
been described and demonstrated. 

A creep crack growth simulation for 304 stainless steel has shown that 
for realistic load levels and corresponding crack speeds the crack-tip field 
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is essentially at a steady-state creep condition. This means that for this 
material the propagating crack-tip field is largely unaffected by the history 
of crack-growth or the history of loading. This feature can greatly reduce 
the analysis requried for predicting creep crack growth behavior in a component 
as can be seen from the following suggested methodology. 

We assume that the crack propagation speed -rr is related to (T.) (i.e., 

01 U X OSS 

- -^) through the poer law suggested by experimental data [13,14] 


^=a[(T,) (11.27) 

dt 1 css 

Next we determine (eg. by steady-state creep finite element analysis) 
(Tj^)^^g as a function of crack length. Because of the assumed steady-state 
crack-tip behavior, this can be accomplished by considering several discrete 
crack lengths and then fitting a curve. No crack growth simulation procedures 
are necessary. Combining (11.27) with this result provides the followng rela- 


tionship between time and crack length 
a(t) 

t = r a (11.28) 

•' ^o 

where a is the initial crack length and t. is the time when crack growth 
initiates. The only unknown quantity in (11.28) is the initiation time t^. 

Vitek [23] has simulated several experiments (compact and double-edge- 
crack specimens) on two CrMoV steels using a dislocation model and has conclu- 
ded that COD correlated well with the initiation of crack growth in these ex- 
periments. If the same conclusion is valid for 304 stainless steel, then one 
can presumably predict t^ based on a transient finite element analysis of the 
initial flawed configuration and a critical value of COD. If initiation occurs 
long after steady-state conditions are reached, it is then reasonable 
to estimate t using the rate of COD obtained from a steady-state finite ele- 
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ment solution. At this time, the validity of (11.28) and of the critical COD 
concept has not been investigated by the authors. 
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APPENDIX A 


Plane Strain ; 

The deformation studied in the example accompanying Part I of this re- 
port is plane strain in the character. Just as for formulations using ordinary 
stresses, a number of the components of the velocity, spin, and stress rate vanish 

if a Cartesian coordiante system is chosen with one axis normal to the plane of 

2 

deformation. We have chosen the x coordinate line to be normal to the plane 
of deformation, so that the velocity, spin, stress rate, and stress are of 
the forms 

= V + V 

13 ^ 31 

m = 0) ^ —3—1 


: *11 rl3 . ;22 

t = t + t e^e3 + t e^ 

^ -31 ^ *13 

+ t e.e. + t e_e_ 

—3—1 — 3r^3 

11 ,13 , 22 

X = T e j^e^ + X e^e 3 + x e^ 

^31 ^33 

+ ^303 + X 0303 

2 

None of the components depends upon x . The velocity is represented on each 
element as 

NQ 

V = E N.q 
1=1 ^ 


The shape functions ^ are described below. Similarly the spin and stress rate 
are represented as 
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0) 


NW . 

= E QW a ; t 
1=1 ' 


NT 

E 

1=1 


qr^S 


We note that this approach requrles minimal specialization In programming 
for the particular case of plane strain. The plane strain condition Is not satis- 
fied a priori; that Is 

5 £22 ^ ° 


for arbitrary 6r. Rather, ^22~^ follows from the stationary condition (a 
component of 7 .1) : 





In using the finite element algorithm the plane strain condition Is only satis- 
fied approximately. In practice a qulltatlve check for satisfaction of the 

22 

plane strain condition can be made by seeing that the stress component x and 
the mean stress are nearly equal. This method for checking '^ 22 ”^ works so long 
as the Inelastic stretching Is proportional to the stress deviator (In the con- 
stitutive equation) . 

SHAPE FUNCTIONS FOR VELOCITY, STRESS RATE, AND SPIN 


Shape Functions for Plane Strain 


x^=X 


3 

X =z 


VELOCITY SHAPE FUNCTIONS 


Ni=N^,.ei+N3,^^^ 


FOUR NODED ELEMENT: 


(l+c?^)(l+nnj 1=1, 2, 3, 4. 

'0 1=5,6, 7, 8. 
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0 


i=l,2,3,4. 





i=5,6,7,8. 


kill, kill, 

= -1, ?2=i. 

n± = - 1 , 12=-!. H3=i, n^=i 

SHAPE FUNCTIONS FOR SPIN: 

QW(1,3,1) = 

QW(3,1,1) = -C^ 

QW(1,3,2)=X C 2 

QW(1,3,3)=Z C 3 

QW(3,1,2) = -X 

QW(3,1,3)=-Z 




31,1-3-1 


The constants were used to improve the 
condition of [H], 


STRESS SHAPE FUNCTIONS: 
QT(1,1,1)=1 
QT(3,1,2)=-1 
QT(2,2,3)=1 
QT(1,3,4)=-1 
QT(3,3,5)=1 
QT(1,1,6)=X 
QT(3,1,6)=-Z 
QT(3,1,7)=-X 
QT(2,2,8)=X 
QT(1,3,9)=-X 
QT(3,3,9)=Z 
QT(3,3,1)=X 


+ 0 + QT22_i£2i2 + '’’^33 . 1-^323 
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QT(1,1,11)=Z 

QT(1,3,12)=-Z 

QT(2,2,13)=Z 

3x3 Gauss Quadrature on this Element. 
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Figure 1. Finite Element Mesh and Boundary Conditions for the 
Problem of Growth of a Void. 
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L=1.00 


L=1.10 


L=1.25 



L=1.35 L=1.50 


Figure 6. Deformation of Cell. 
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The 102 element mesh (331 nodes; 642 d.o.f.) 
The 300 element mesh (941 nodes; 1840 d.o.f) 



Figure 13, Summary of Geometry, Loading, Material Properties and 
Finite Element Meshes for the Compact Specimen. 
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I 


Properties Representative of 
304 Stainless Steel at 650®C 


s,s 




E = 1.5 X 10® MPa 
V =0.3 

y = 4 X lO"'® MPa‘>hr 
n = 7 


Figure 15. 


Sununary of Geometry, Loading, Material Properties and Finite Element Meshes for 
the Compact Specimen. 
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Figure 16. 



Illustration of Mesh Shifting/Remeshing Procedure for 
Simulation of Crack Grov/th. 
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Figure. 23. Values of (f) and C* During 

Earlier Phases of Creep Crack 
Growth . 




Figure 24. Values of (T) and C* During Latter Phases of Creep 
Crack Growth. 


365 


